package com.example.analysis.utils;

import com.example.analysis.exception.BusinessException;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.LUDecomposition;
import java.util.Arrays;

public class ComCalcUtils {
    /**
     * 计算均值
     *
     * @param data 数据
     * @return 均值
     */
    public static double calcMean(double[] data) {
        if (data.length == 0)
            return 0;
        double mean = 0.0;
        for (Double item : data) {
            mean += item;
        }
        mean /= data.length;
        return mean;
    }

    /**
     * 计算标准差
     *
     * @param data 数据
     * @return 标准差
     */
    public static double calcStd(double[] data, double mean) {
        if (data.length == 0)
            return 0;
        double sumSq = 0;
        for (double x : data) sumSq += Math.pow(x - mean, 2);
        return Math.sqrt(sumSq / (data.length - 1));
    }

    // 计算中位数
    public static double calcMedian(double[] data) {
        Arrays.sort(data);
        int n = data.length;
        if (n % 2 == 0) {
            return (data[n / 2 - 1] + data[n / 2]) / 2.0;
        } else {
            return data[n / 2];
        }
    }

    /**
     * 计算方差
     */
    public static double calcVariance(double[] data) {
        double avg = calcMean(data);
        return Arrays.stream(data).map(x -> Math.pow(x - avg, 2)).average().orElse(0);
    }

    // z-score标准化
    public static double[][] zScore(double[][] data) {
        int n = data.length, m = data[0].length;
        double[][] normalizedData = new double[n][m];
        for (int i = 0; i < n; i++) {
            // 计算均值
            double mean = calcMean(data[i]);
            // 计算标准差
            double std = calcStd(data[i], mean);
            for (int j = 0; j < m; j++) {
                // 进行 Z-score 标准化
                normalizedData[i][j] = (data[i][j] - mean) / std;
            }
        }
        return normalizedData;
    }

    // 极值法，正向化处理数据/归一化（消除量纲）/标准化（消除数量级）
    public static void minmax(double[] data) {
        double max = Arrays.stream(data).max().orElse(0), min = Arrays.stream(data).min().orElse(0), denominator = max - min;
        for (int i = 0; i < data.length; i++) {
            data[i] = (data[i] - min) / denominator;
        }
    }

    // 负向化处理数据
    public static void n_minmax(double[] data) {
        double max = Arrays.stream(data).max().orElse(0), min = Arrays.stream(data).min().orElse(0), denominator = max - min;
        for (int i = 0; i < data.length; i++) {
            data[i] = (max - data[i]) / denominator;
        }
    }

    // 适度化
    public static void moderate(double[] data) {
        double max = Arrays.stream(data).max().orElse(0), min = Arrays.stream(data).min().orElse(0), denominator = max - min, mean = calcMean(data);
        for (int i = 0; i < data.length; i++) {
            data[i] = Math.abs(data[i] - mean) * -1;
            data[i] = (data[i] - min) / denominator;
        }
    }

    // 计算Pearson系数
    public static double calcPearson(double[] x, double[] y) {
        // 均值
        double xMean = ComCalcUtils.calcMean(x);
        double yMean = ComCalcUtils.calcMean(y);
        // 计算Pearson系数
        double numerator = 0, denominator, a = 0, b = 0;// 分子，分母
        for (int i = 0; i < x.length; i++) {
            numerator += (x[i] - xMean) * (y[i] - yMean);
            a += Math.pow((x[i] - xMean), 2);
            b += Math.pow((y[i] - yMean), 2);
        }
        denominator = Math.sqrt(a) * Math.sqrt(b);
        if (denominator == 0)
            return 0;
        return numerator / denominator;
    }

    // ---计算相关系数矩阵
    public static double[][] calculateCorrelationMatrix(double[][] data) {
        int p = data.length;    // 现在 p 是变量的数量
        double[][] correlationMatrix = new double[p][p];
        for (int i = 0; i < p; i++) {
            for (int j = 0; j < p; j++) {
                correlationMatrix[i][j] = calcPearson(data[i], data[j]);
            }
        }
        return correlationMatrix;
    }

    // 矩阵求逆
    public static double[][] inverseMatrix(double[][] matrix) {
        int n = matrix.length;
        double[][] augmented = new double[n][2 * n];

        // 构建增广矩阵 [A | I]
        for (int i = 0; i < n; i++) {
            System.arraycopy(matrix[i], 0, augmented[i], 0, n);
            augmented[i][i + n] = 1;
        }
        // 高斯 - 约旦消元法求逆
        for (int i = 0; i < n; i++) {
            double pivot = augmented[i][i];
            for (int j = 0; j < 2 * n; j++) {
                augmented[i][j] /= pivot;
            }

            for (int k = 0; k < n; k++) {
                if (k != i) {
                    double factor = augmented[k][i];
                    for (int j = 0; j < 2 * n; j++) {
                        augmented[k][j] -= factor * augmented[i][j];
                    }
                }
            }
        }
        // 提取逆矩阵
        double[][] inverse = new double[n][n];
        for (int i = 0; i < n; i++) {
            System.arraycopy(augmented[i], n, inverse[i], 0, n);
        }
        return inverse;
    }

    // 计算偏相关系数矩阵
    public static double[][] calculatePartialCorrelationMatrix(double[][] correlationMatrix) throws BusinessException {
        // 1.是否可逆
        //  1.1 转成矩阵
        Array2DRowRealMatrix corrMatrix = new Array2DRowRealMatrix(correlationMatrix);
        //  1.2 计算相关矩阵的行列式（需处理奇异矩阵）
        LUDecomposition luDecomp = new LUDecomposition(corrMatrix);
        double det = luDecomp.getDeterminant();
        if (det <= 0) {
             throw new BusinessException("相关矩阵不可逆，无法进行KMO检验");
        }

        int p = correlationMatrix.length; // 题项数量
        double[][] inverse = inverseMatrix(correlationMatrix);
        double[][] partialCorrelationMatrix = new double[p][p];

        for (int i = 0; i < p; i++) {
            for (int j = 0; j < p; j++) {
                if (i != j) {
                    partialCorrelationMatrix[i][j] = -inverse[i][j] / Math.sqrt(inverse[i][i] * inverse[j][j]);
                }
            }
        }
        return partialCorrelationMatrix;
    }

    // ---计算协方差矩阵
    public static double[][] calculateCovarianceMatrix(double[][] data) {
        int n = data[0].length; // 样本数量
        int p = data.length;    // 题项数量
        double[][] covarianceMatrix = new double[p][p];
        for (int i = 0; i < p; i++) {
            double meanI = 0;
            for (int k = 0; k < n; k++) {
                meanI += data[i][k];
            }
            meanI /= n;

            for (int j = 0; j < p; j++) {
                double meanJ = 0;
                for (int k = 0; k < n; k++) {
                    meanJ += data[j][k];
                }
                meanJ /= n;

                double sum = 0;
                for (int k = 0; k < n; k++) {
                    sum += (data[i][k] - meanI) * (data[j][k] - meanJ);
                }
                covarianceMatrix[i][j] = sum / (n - 1);
            }
        }
        return covarianceMatrix;
    }
}
